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LONG-TERM  GOALS 

Our  goal  is  to  provide  a  complete  simulation  code  that  will  represent  and  predict  the  sediment  transport 
and  bed  features  on  the  continental  shelf  at  user-specified  resolution  by  using  state-of-the-art  algorithms 
for  the  physics  and  numerics  of  the  simulation  code. 

OBJECTIVES 

Our  primary  objective  is  to  simulate  the  ripple  climate  on  the  bed  of  the  inner  shelf  at  depths  on  the  order 
of  20  m  and  over  domains  ranging  from  centimeters  to  kilometers  in  support  of  the  SAX  experiments 
and  analyses.  Our  secondary  objective  is  to  use  simulation  to  better  understand  the  physics  of  ripple 
formation  and  sediment  transport  in  this  environment.  A  third  objective  is  to  create  a  powerful  and  fast 
parallel  numerical  code  system  by  the  synthesis  of  a  number  of  our  computational  tools. 

APPROACH 

The  Critical  Benthic  Environmental  Processes  and  Modeling  at  SAX04  (aka  Ripples)  Departmental  Re¬ 
search  Initiative  of  ONR  has  a  number  of  objectives.  Among  them  our  work  is  addressing  the  following: 

•  Model  ripple  morphology  and  its  gradients  on  scales  ranging  up  to  kilometers.  Model  the  effects 
of  bioturbation  on  the  ripple  evolution. 

•  Understand  the  response  of  ripples  to  changes  in  wave  and  wave-current  forcing. 

Based  on  our  experience  in  oceanic  and  atmospheric  field-experiment  programs,  we  also  understand  the 
importance  of  the  close  and  timely  collaboration  of  all  of  the  investigators  in  the  modeling  and  field 
components  of  the  program.  For  the  DRI  we  are  aiming  to  be  able  to  make  predictive  simulations  of  the 
bed  evolution  for  the  second  field  season  in  FY06.  We  are  using  the  conditions  from  SAX04  as  the  basis 
for  hindcasts  and  code  test  and  verification  using  the  synthesis  of  codes  described  below. 


1 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

30  SEP  2005 


2.  REPORT  TYPE 


3.  DATES  COVERED 

00-00-2005  to  00-00-2005 


4.  TITLE  AND  SUBTITLE 

Simulation  of  Benthic  Ripples  and  Transport  Processes  for  SAX 


6.  AUTHOR(S) 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 
5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES)  8.  PERFORMING  ORGANIZATION 

Stanford  University, Stanford, CA,  94305  report  number 

9.  SPONSORING/MONITORING  AGENCY  NAME(S )  AND  ADDRESS(ES )  10.  SPONSOR/MONITOR' S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

code  1  only 

14.  ABSTRACT 

Our  primary  objective  is  to  simulate  the  ripple  climate  on  the  bed  of  the  inner  shelf  at  depths  on  the  order 
of  20  m  and  over  domains  ranging  from  centimeters  to  kilometers  in  support  of  the  SAX  experiments  and 
analyses.  Our  secondary  objective  is  to  use  simulation  to  better  understand  the  physics  of  ripple  formation 
and  sediment  transport  in  this  environment.  A  third  objective  is  to  create  a  powerful  and  fast  parallel 
numerical  code  system  by  the  synthesis  of  a  number  of  our  computational  tools. 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

8 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


Field-scale  solver  We  are  using  the  parallel  unstructured  grid  code,  known  as  SUNTANS  (Stanford 
Unstructured  Nonhydro  static  Terrain-following  Adaptive  Navier-Stokes  Simulator),  which  solves  the 
nonhydro  static  Navier-Stokes  equations  under  the  Boussinesq  approximation  with  the  large-eddy  for¬ 
mulation  for  the  resolved  motions.  The  code  has  been  developed  by  Fringer  et  al.  (2005)  in  order  to 
determine  the  behavior  of  the  nonhydrostatic  internal  wave  spectrum  in  Monterey  Bay  (Fringer  et  al., 
2004).  SUNTANS  is  a  high-resolution  simulation  tool  for  coastal,  estuarine,  and  river  simulations,  and 
incorporates  a  host  of  processes,  from  the  free-surface  to  the  bed,  that  play  an  important  role  in  driving 
smaller-scale  motions,  including  the  nonhydrostatic  pressure,  wind  stress,  surface  and  internal  tides,  and 
parameterizations  of  small-scale  internal  as  well  as  surface  breaking  physics.  Most  importantly,  because 
it  is  unstructured,  and  adaptive  in  the  long-term,  we  have  the  ability  to  employ  higher  resolution  around 
the  site  of  interest  for  the  SAX  field  experiments  where  we  embed  the  smaller-scale  simulation  codes 
that  we  describe  below. 

Parallel  large-eddy  simulation  code  The  SUNTANS  code  is  being  used  to  determine  the  field-scale 
motions  that  drive  the  smaller-scale  flow,  which  is  computed  with  a  nonhydro  static  parallel  large-eddy 
simulation  Navier-Stokes  solver.  This  solver  employs  the  code  of  Cui  and  Street  (2001;  2004)  (hence¬ 
forth  referred  to  as  PCUI)  which  is  an  adaptation  of  the  large-eddy  simulation  solver  of  Zang  et  al. 
(1994)  to  simulate  laboratory-scale  rotating  stratified  flows  using  MPI,  the  message  passing  interface, 
on  parallel  computers.  The  formulation  of  Zang  et  al.  employs  a  fully  non-orthogonal  curvilinear, 
boundary-following  grid  as  well  as  a  dynamic-mixed  sub  filter  scale  model  (Zang  et  al.,  1993)  to  com¬ 
pute  the  subfilter  scale  terms  that  arise  from  volume  filtering.  The  equations  are  discretized  in  time  with  a 
semi-implicit,  second-order  accurate  method  and  in  space  with  second-order  differencing  and  advanced 
with  a  fractional  step/projection  method.  The  method  has  been  applied  to  a  host  of  simulations,  including 
turbulent  stratified  flow  over  a  wavy  bed  (Calhoun  et  al.,  1999),  interfacial  breaking  waves  (Fringer  & 
Street,  2003),  and  the  study  of  a  round  jet  in  crossflow  (Yuan  et  al.,  1999)  using  a  total  of  20  overlapping 
grids.  Along  with  having  been  verified  extensively,  the  PCUI  code  is  highly  efficient,  as  it  was  41% 
faster  than  the  next  comparable  code  on  the  NASA  IBM  SP-2  (Bergeron,  1998). 

Sediment  transport  For  sediment  transport,  we  are  employing  algorithms  from  the  code  of  Zedler  and 
Street  (2001;  2002;  2005),  who  extended  the  formulation  of  Zang  et  al.  (1994)  and  added  a  module 
which  solves  the  advection  diffusion  equation  (with  settling  term)  for  the  sediment.  Under  this  work 
Zedler  and  Street  also  created  a  bed-load  module  and  a  bed  erosion/deposition  module  (Zedler,  personal 
communication,  2005)  to  go  with  the  module  for  suspended  sediment  transport.  She  provided  us  with  the 
equations  in  the  useful  terrain-following  coordinates  of  her  original  code  so  that  the  resolution  of  shear 
stress  direction,  for  example,  is  clear. 

Motion  of  the  bed  may  be  computed  using  the  immersed  boundary  method  (IBM),  which  allows  the 
tracking  of  complex  boundaries  on  Cartesian  grids  (Tseng  &  Ferziger,  2003),  or  by  the  moving  bound¬ 
ary  method  of  Hodges  and  Street  (1999)  and  Fringer  (2003)  which  preserves  the  near-bed  terrain  fol¬ 
lowing  coordinates.  The  IBM  method  effectively  imposes  forcing  functions  on  the  right  hand  side  of 
the  momentum  equations  that  enforce  boundary  conditions  at  complex  water-ripple  boundaries;  under 
the  related  project  cited  below,  Zedler  has  created  a  preliminary  implementation  of  an  IBM  for  flow  and 
sediment  transport  at  a  bed. 

The  Zedler  code  was  a  serial  code  and  large  problems  ran  very  slowly;  thus,  the  elements  of  the  Zedler 
and  Street  code  are  currently  being  implemented  into  the  PCUI  code  to  run  in  parallel  mode.  This  code 
is  referred  to  as  PCUUSed.  Both  options  for  moving  the  boundary  are  being  kept  open  at  this  time. 
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Strategy:  We  are  molding  the  above  tools  into  an  integrated  system  wherein: 


•  An  enhanced  SUNTANS  drives  the  sediment  transport  code  PCULSed  to  simulate  local  ripple 
climate  wherever  that  information  is  needed.  This  need  not  be  done  in  real  time,  but  can  be  done 
rapidly  because  the  codes  are  parallelized. 

•  The  parallelized,  IBM  sediment  transport  code  (PCUI  Sed)  simulates  the  local  bed  evolution.  Our 
plan  is  to  run  the  nested  PCULSed  in  a  domain  of  roughly  10  m  by  10  m  with  5  cm  horizontal 
resolution.  Results  from  this  domain  could  in  principle  drive  a  50  cm  by  50  cm  domain  with  25 
mm  horizontal  resolution.  However,  based  on  SAX99  information,  we  anticipate  a  focus  on  ripple 
lengths  ranging  from  about  0.5  m  down  to  7  cm. 


WORK  COMPLETED 

This  grant  began  in  January  2005.  At  this  time  we  are  in  the  process  of  implementing  the  immersed 
boundary  method  into  SUNTANS  to  create  a  terrain-following  boundary.  We  report  on  this  implementa¬ 
tion  below.  Our  next  steps  involve  modification  of  the  PCUI  parallel  code  to  create  PCULSed  using  the 
tools  created  by  Zedler  under  the  related  project  cited  below. 

RESULTS 

The  immersed  boundary  method  employs  the  use  of  Cartesian  grids  in  complex  geometries  by  adjusting 
ghost  cell  velocity  and  pressure  values  in  such  a  way  that  boundary  conditions  are  satisfied  on  the  real, 
or  immersed,  physical  boundaries,  rather  than  at  the  faces  of  the  Cartesian  grid.  The  ghost-cell  velocities 
and  pressure  field  are  extrapolated  from  the  computational  cells  inside  the  immersed  boundary  while 
satisfying  boundary  conditions  at  that  boundary.  On  staggered  grids,  which  are  employed  in  SUNTANS, 
the  extrapolation  procedure  depends  on  whether  the  immersed  boundary  lies  between  vertical  or  hori¬ 
zontal  faces,  and  it  is  quite  different  for  the  two.  For  horizontal  faces,  the  computation  in  the  absence 
of  the  immersed  boundary  assumes  that  the  bottom  boundary  coincides  with  the  face  of  a  Cartesian  grid 
cell,  as  depicted  in  Figure  1(a).  Inclusion  of  the  immersed  boundary  effectively  alters  the  position  of  the 
no-slip  boundary  condition  and  as  a  result  the  value  of  the  ghost-cell  velocity  is  altered  to  account  for 
this,  as  depicted  in  Figure  1(b).  For  vertical  faces,  the  horizontal  velocity  vectors  at  the  vertical  faces 
of  the  Cartesian  boundary  are  zero  in  the  absence  of  the  immersed  boundary,  while  its  inclusion  adds 
nonzero  velocities  to  these  faces.  This  in  turn  alters  the  horizontal  and  vertical  velocity  field  to  account 
for  the  complex  geometry  more  accurately. 

The  actual  implementation  of  the  immersed  boundary  method  is  referred  to  in  the  literature  as  ghost 
velocity  reconstruction  (Tseng  &  Ferziger,  2003).  Using  the  procedure  of  Tseng  and  Ferziger,  the  ghost 
values  beneath  the  bottom  boundary  are  computed  through  linear  extrapolation,  which  is  advantageous 
on  the  unstructured,  z-level  grids  in  SUNTANS  (see  Figure  4  for  a  schemetic  of  the  z-level  prism  cells) 
because  it  only  requires  velocities  on  the  upper  and  lower  sides  of  the  immersed  boundary  rather  than 
at  other  faces  of  the  unstructured  grid,  which  would  require  complicated  and  impractical  interpolation 
procedures.  Using  linear  extrapolation  in  the  vertical  has  the  added  advantage  of  guaranteeing  that  the 
ghost  velocity  is  normal  to  the  vertical  face  where  it  is  located. 

Following  the  work  of  Kim  et  al.  (2001),  the  construction  method  at  bottom  boundaries  is  illustrated 
in  Figure  3.  Here,  U  represents  the  component  of  the  velocity  normal  to  vertical  faces,  while  d0, 
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Figure  1:  Depiction  of  the  immersed  boundary  method  as  applied  to  horizontal  faces  showing  its 
effect  on  the  ghost  velocities  (in  bold)  without  (a)  and  with  (b)  the  immersed  boundary. 


Figure  2:  Depiction  of  the  immersed  boundary  method  as  applied  to  vertical  faces  showing  its  effect 
on  the  ghost  velocities  (in  bold)  without  (a)  and  with  (b)  the  immersed  boundary.  The  shaded  region 
shows  that  the  immersed  boundary  also  indirectly  affects  the  vertical  velocity  field  in  the  cells 

containing  the  immersed  boundary. 


and  d2  represent  the  distance  between  the  immersed  boundary  (IB)  and  the  ghost  point,  the  center  of 
the  bottommost  cell  inside  fluid  domain,  and  the  center  of  the  cell  right  above  the  bottommost  cell, 
respectively.  In  order  to  ensure  stability,  the  ghost  velocity  must  be  bounded  by  the  minimum  and 
maximum  values  of  the  velocity  field  within  the  immersed  boundary.  This  is  ensured  if  the  interpolation 
is  altered  when  the  distance  between  the  ghost  velocity  and  the  immersed  boundary  exceeds  d\,  such 
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that,  if  Ug  is  the  ghost  velocity  and  Uk  and  Uk-i  are  the  velocities  above  the  immersed  boundary,  then 


Un 


—  ^h^k  0<do<di, 

(d*-d0)uk+(do-d1)uk-i  otherwiSe . 

a  2  — «i 


(a) 


(b) 


Figure  3 :  Depiction  of  the  immersed  boundary  method  of  obtaining  the  ghost  velocity  Ugfor  the  case 
when  (a)  the  ghost  velocity  is  closer  to  the  immersed  boundary  than  the  computational  velocity  and 
(b)  the  ghost  velocity  is  farther  away  from  the  immersed  boundary  than  the  computational  velocity. 

Because  SUNTANS  grids  are  unstructured  in  plan,  the  formulation  is  much  more  complex  for  the  vertical 
faces,  since  boundary  conditions  must  be  satisfied  at  piecewise  planar  immersed  boundaries  that  intersect 
prism  cells,  as  depicted  in  Figure  4.  The  method  is  similar  to  that  depicted  in  Figure  3  except  the 
interpolation  is  much  more  complex  because  of  the  unstructured  nature  of  the  horizontal  interpolation. 

While  the  velocity  field  must  satisfy  Dirichlet  conditions  at  the  immersed  boundaries,  the  pressure  field 
is  updated  by  satisfying  a  Neumann  condition,  since  there  must  be  a  zero  pressure  gradient  normal  to 
immersed  boundaries.  This  is  naturally  satisfied  for  the  pressure  field  because  it  is  specified  at  the  centers 
of  the  cells  rather  than  at  the  faces,  which  complicates  the  immersed  boundary  method  for  the  velocity 
field.  To  compute  the  value  of  the  pressure  field  in  the  ghost  cells,  we  assume  the  immersed  boundary 
that  intersects  an  unstructured  grid  cell  is  given  by  the  piecewise-planar  function  (given  by  the  red  lines 
in  Figure  4) 

ax  +  by  +  cz  —  d ,  (1) 

such  that  \/ a2  +  b'2  +  c?  =  1.  Following  Tseng  and  Ferziger  (2003),  if  the  value  of  the  pressure  in  the 
vicinity  of  the  boundary  is  assumed  to  be  of  the  form 

p(x,  y,  z)  =  wxx  +  w2y  +  w3z  +  w4  ,  (2) 

then  it  follows  that  the  value  of  the  pressure  at  the  ghost  point  in  Figure  5  is 

Pg  =  WiXg  +  W2yg  +  W3Zg  +  W*  ,  (3) 
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while  the  pressure  gradient  is  VP  =  w±ex  +  w 2ey  +  w 3ez,  where  ex,y,z  are  the  Cartesian  unit 
vectors.  If  represents  the  normal  to  the  plane  defined  in  equation  (1),  then  the  component  of  the 
pressure  gradient  in  the  direction  of  nb  is  given  by 

dp 

—  —  VP  •  rib  —  aw i  +  bw 2  +  cw 3  . 

On  b 

Using  this  equation,  along  with  three  constraints  that  require  the  pressure  at  points  1,  2,  and  3  in  Figure 
5  to  satisfy  equation  (2),  one  arrives  at  the  governing  equations  for  the  weights  in  equation  (3)  as 
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Here,  pi,  P2,  and  p3  are  the  pressure  in  the  computational  cells  inside  the  immersed  boundary  which  are 
located  at  ( x ,  y,  z)  1,  ( x ,  y,  z)2,  and  ( x ,  y,  z) 3,  respectively,  as  in  Figure  5,  and  ( dp/ dn)b  =  0  is 
the  desired  pressure  gradient  normal  to  the  immersed  boundary. 


(b) 


Figure  4:  Depiction  of  how  the  bottom  topography  (in  gray )  intersects  the  triangular  prisms  in 
piecewise-planar  faces  (red),  (a)  Three-dimensional  view,  (b)  top  view,  showing  the  velocity  vector 
defined  at  the  cell  center  u0  and  the  velocity  interpolated  to  the  immersed  boundary  edge  ub. 


RELATED  PROJECTS 

This  project  is  linked  to  the  just  completed  ONR  Grant  N000 14-00- 1-0440  “Large  Eddy  Simulation  of 
Sediment  Transport  in  the  Presence  of  Surface  Gravity  Waves,  Currents  and  Complex  Bedforms”  from 
which  we  are  adapting  algorithms  created  by  Dr.  Zedler  for  sediment  transport  and  bed  movement. 
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Figure  5:  Vertical  slice  depicting  how  the  normal  pressure  gradient  condition  is  satisfied  at  the 
immersed  boundary  using  the  pressure  at  three  computational  cells.  The  red  line  represents  the 
intersection  of  the  red  plane  in  Figure  4  with  the  vertical  plane  represented  here. 
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